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ABSTRACT 

We examine the temperature structure of the IGM due to the passage of individual 
ionisation fronts using a radiative transfer (RT) code coupled to a particle-mesh (PM) 
iV-body code. Multiple simulations were performed with different spectra of ionising 
radiation: a power law (oc i/~ ' 5 ), miniquasar, starburst, and a time- varying spectrum 
that evolves from a starburst spectrum to a power law. The RT is sufficiently resolved 
in time and space to correctly model both the ionisation state and the temperature 
across the ionisation front. We find the post-ionisation temperature of the reionised 
intergalactic medium (IGM) is sensitive to the spectrum of the source of ionising 
radiation, which may be used to place strong constraints on the nature of the sources 
of reionisation. Radiative transfer effects also produce large fluctuations in the Hen 
to Hi number density ratio r/. The spread in values is smaller than measured, except 
for the time- varying spectrum. For this case, the spread evolves as the spectral nature 
of the ionising background changes. Large values for r\ are found in partially ionised 
Hen as the power-law spectrum begins to dominate the starburst, suggesting that the 
large 77 values measured may be indicating the onset of the Hen reionisation epoch. 

Key words: radiative transfer - cosmology: diffuse radiation - cosmology: large-scale 
structure of Universe - methods: numerical - methods: N-body simulations 



^ ; 1 INTRODUCTION 



ies or st ar clusters) |Choudhurv fc Ferraral 120061) or mini- 



The process by which the Universe was reionised is one of 
the premier unsolved questions in cosmology. Measurements 
of the Lya optical depth of the Intergalactic Medium (IGM) 
along Quasi-Stellar Object (QSO) lin es-of-sight require th e 
IGM to have been reionised by z>6 (|Becker fc et alll200ll ). 
consistent with recent measurements of the Cosmic Mi- 
crowave Background (CMB) by the Wilkinson Microwave 
Anisotropy Probe (WMAP), although a higher redshift of 
z ~ 11 is preferred , with a 2-a upper limit of z < 17 
l|Spergel fc et alll2007l ). The sources of the reionisation are 
currently unknown. The most recent estimates of the num- 
bers of high redshift QSOs suggest QSOs are too few to 
have ionised the Hi prior to z « 4, without an upturn 
in th e QSO luminosity function at the faint end l|Meiksinl 
120051 ). While an adequate supply of ionising photons is 
likely produced within young, star-forming galaxies, it has 
yet to be conclusively demonstrated that the ionising pho- 
tons are able to escape in suf ficient numbers to meet the 
requirements for re ionisation l|Fernandez-Soto et al.l 120031 ; 
iMalkan et ail 12003 V Other, more speculative, possibilities 
include pockets of Population III stars (e.g., in young galax- 
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quasars l|Madau et al.ll2004l ) 



The post-ionisation IGM temperature may provide a 
clue to the nature of the reionisation sources. Low density 
regions retain a memory of the post-ionisatio n temperature 
l|Meiksinlll993 : lMiralda-Escude fc Reeslll994l ). Evidence for 
temperatures in excess of the optically thin predictions at 
high redshifts [z > 3) is provided by t he widths of the low 
column density Lya forest absorbers (M eiksin et al.ll200ll ) 
which are expected to reside in the underdense regions com- 
prisin g most of the volume of the Universe (|Zhang et alj 
|199S() . 

In the past few years, several groups have implemented 
numerical radiative transfer schemes to solve for the reion- 
isation of the IGM (lAbel et al.|[l99i iGnedin fc Abelll200ll; 
Nakamoto et alj l200ll: iRazoumov et al.1 120021: ICiardi et all 



20031 : llliev et all 120061 : IWhalen fc Norman! 120061 ) . Most of 



the modelling of IGM reionisation has focused on the prop- 
agation of ionisation fronts (I-fronts), with less emphasis 
given to t he accurate com putation of post-ionisation tem- 
peratures (|lliev et al.ll2006l ). as the requi red calculations are 
much more computationa lly demanding l|Bolton et al.|[20o3 : 
iBolton fc Haeh nclt 2007). Crucial to accurate temperatures 
is the resolution of the ionisation fronts, particularly for hard 
spectrum sources like young galaxies or Active Galactic Nu- 
clei (AGN) . Early estimates of the post-ionisation IGM tem- 
perature were made assuming optically thin gas. Accounting 
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for radiative transfer within the ionisation fronts produces 
much higher temperatures. 

In this paper, we focus on the thermal properties of 
the post-ionised IGM. We investigate the role played by 
the spectral shape of the dominant ionising sources on the 
temperature of the post-ionised IGM, considering both soft 
sources like galaxies and hard like QSOs. Rather than per- 
forming a full reionisation simulation with multiple sources, 
which is still beyond the current resources of numerical com- 
putations, we limit the computations to a single ionisation 
front sweeping across the box, as would occur prior to the 
time of the complete overlap of ionisation fronts. Although 
this approximation will not allow us to compute precise val- 
ues for the ionisation fractions of hydrogen and helium for 
comparison with measurements, as these levels will be re- 
set by the subsequent total photoionisation rates after I- 
fronts overlap, it does permit a quantitative evaluation of 
the temperature structure after complete reionisation. This 
is because, once ionised, the temperature is insensiti ve to the 
total intensity or shape of the ionising background l|Meiksinl 
ll994l ; lHui fc Gnedinlfi997T l. which will readjust both due to 
the overlapping of I-fronts and the evolution of the sources 
and photoelectric optical depth of the IGM. 

The simulations are performed with a radiative transfer 
code coupled to an iV-body code to model the evolution of 
the IGM. In ij2] the algorithm and the numerical code are 
described. The details of the simulation volume and sources 
are given in SJ3] Results of the simulations are provided in 
£|4] and discussed in fj5] 



2 METHOD 

The simulation code, PMRTjnpi, is the merger of a La- 
grangian particle-mesh (PM) code (I Meiksin et al.i 1999), and 
a grid-based radiative transfer (RT) code which are mod- 
ularly independent. We perform the RT on the evolving 
gas density field, as computed by the PM code. Full hy- 
drodynamical simulations have shown that the gas den- 
sity closely traces the dark matter density down to th e 
Jeans scale (~ 100 kpc) (|Cen et al.ll 19941 ; IZhang et al.ll 19981 ). 
Statistical comparisons between the resulting Lycv forest 
iroperties show good agreement at the 10 per cent level 
Meiksin fc Whitell200ll ; lviel et al.ll2002l ). 

The gas density is taken as a constant fraction (£l),/Q m ) 
of the total density. This is a fair appro ximation since 
baryons trace dark matter (|Zhang et al.lll998r ) and in simu- 
lations with gas, the gas to dark matter ratio doesn't vary 
more than a facto r of 30 percent over the over density range 
S = 0.1 to 1000 l|Tittlev fc Couchmanl koool) . The 1-fronts 
propagate much faster than the sound speed, so that the 
pressure-response of the gas has only a small effect on the 
reionisation. The PM code evolves the density field assuming 
all mass is collisionless, interacting only gravitationally. 

The density field is determined from the particle distri- 
bution by griding the particles on to a mesh. To avoid low- 
count artefacts in low-density regions while not sacrificing 
information in dense regions, the density field is separated 
into two fields, p(r ) = pio(r) + phi(r), and the low-density 
field, pi a (r), is convolved with a Gaussian of radius two grid 
cells. A threshold cell count (N t h res h = 5 for these simu- 
lations) discriminates pio(jc) from phi(r). In dense regions, 



Table 2. Ionisation cross-section parameters used in Eq. \E\ 
or[m 2 ] vt [H z ] P s 



ct H i 6.30 x icr 22 
cr HcI 7.83 x 1CT 22 

(THell 1-58 X 1CT 22 



3.282 x 10 15 
5.933 x 10 15 
1.313 x 10 16 



1.34 2.99 
1.66 2.05 
1.34 2.99 



pio(r) = Nthresh and p w (r) = p(r) - N thresh , while in 
low-density regions, pi {r) = p(r). The smoothed field is 
then p (v ) = pio(jc ) ® p(r ) + phi(jc ) where g is the Gaussian 
smoothing kernel. 

2.1 Radiative transfer 

The RT code uses a probabilistic me thod which i s bas ed on 
the photon-conserving alg orithm of | Abel et al.l l| 19991 '). ex- 
tended to include helium bv lBolton et al.l (|2004f ) who applied 
the RT algorithm to a density field frozen in the comoving 
frame. For convenience, we provide details of the RT code 
here. 

The rates of change of the ionisation-state populations 
due to ionisations and recombinations are given by: 

khi = n e nHuctHu — miiTm (1) 
riHii = —fiui 

WHcI = IlelHellClHcII — ^Hc^Hel 

«HcII = ~ nuci — nneiii 

fbHelH = ^HoIiFhcII — ^e^HcIIlQ!HeIII 

where a; is the recombination coefficient from species i and 
Ti is the photoionsation rate. Number-density changes, hi, 
due to cosmological evolution are accounted for by carrying 
only the ionisation fractions (/hi = nui/nn and similar for 
helium fraction /hoi, /hoii, and /Hem) between PM itera- 
tions. 

The forms for the recombination coefficients ami(T) 
and QHe m{T) are d erived f r om th e generalised hydrogenic 
case A form given in Seat onl (1 19591) ForaHe 11 (T) , the radia- 
tive term is from IVerner fc Ferlandl ill99rj) while the dielec- 
tronic term is adopted from Aldrovandi fc Pecmignotl ll 19731 ). 
The recombination coefficients are provided in Table [T] 

The photoionisation rate per particle, T, for each species 
is dependent on the local mean intensity of radiation, J„, 
and the ionisation cross section, a, for the species by 



47T 



hp 



(2) 



where ut is the threshold frequency for ionisation, which 
differs for each species. In the plane wave approximation 
used here, the local mean intensity is simply the radiation 
field incident on the volume attenuated by the cumulative 
optical depth, 



47rJ„ = 



L v 
AnR 2 ' 



(•3) 



where L v is the luminosity of the source per frequency, 7?o 
is the distance of the source from the volume, and t„ is 
the cumulative optical depth which, at distance R from the 
source, is: 



T V (R) — CTHiAhi + OHoiAhcI + OHcIiTVhcII 



(4) 



Reionisation scenarios and the temperature of the IGM 3 



Table 1. Recombination («; [m 3 s x ]) and recombination cooling (/3; [J m 3 s 1 ]) coefficients. 



"II 



„ 2.065 x lQ-^T- 1 /? (6.414 - \ In T + 8.68 x lCr 3 T 1 /3^ 



OHeii 3.294 X 10" 



T ' 
15.54 , 



1 + 



1, -i 0.309 

12 



1 + 



3.676 X10 7 . 



8.260 x 10- 17 T- 1 /2 (V.107 - | In T + 5.47 x 10- 3 tV3^ 



,, , -9.4x10" 

+ 1.9 x lO" 9 1 + 0.3e T 



-T- h 



0Hu 2.851 x lO-^T 1 ^ (5.914 - i In T + 0.01184T 1 /3 

,„ „ „„ / -9.4x10" ~\ -4.7xl0 a 

/3 Ho ii 1-55 x io- 39 T 3647 + 1.24 x 10~ 26 ( 1 + 0.3e T ) e r T~ h 

/3 H oiii 1-140 x 10- 39 TV2 (6.607 - | In T + 7.459 x IO^tVs J 



where the cumulative column depth from the edge of the 
box to R is 

Ni(R) = / ni (r)dr. (5) 

The ionisatio n cross sect i ons a re approximated using the 
form given in lOsterbrockl (|l989l ). 



L = Lhii + LucH + £hoih + L e H + Lc ■ 



(12) 



+ (1-/?) 



(6) 



The parameters for the various species are listed in Table [2] 
The reionisation is further restricted by imposing a check 
on the position of the light front to ensure gas is not ionised 
too early. Neglectin g this effect can lead to an over-estimate 
of the temperature ( Madau et al.|[l997l) . 

The equations governing the thermal state of the gas 
are used in their entropic form. The entropy is parametrized 
by the function 
P 

S = —, (7) 
pi' w 

where P is the pressure, p is the gas density, and 7 = 5/3 
is the adiabatic index for a monatomic gas. It follows from 
Eq. [7] that 

S=(7-1)P~ 7 (G-L), (8) 

where G and L are the thermal gain and loss functions per 
volume, and 



T 



k 



-Sp- 



where fj, is the mean molecular weight of the gas, m u is an 
atomic mass unit, and k is the Boltzmann constant. 

Heating is provided by the excess energy above the ion- 
isation threshold Kvt of the ionising photons. For a single 
species of density n, 



G = 4nn 



(y — vt)<Iv. 



(10) 



Since all species contribute, the total heating rate is 

G = Ghi + Ghci + Ghoii- (11) 

Cooling is provided by recombinations, collisional ex- 
citation of the excited levels in neutral hydrogen, and in- 
verse Compton scattering off cosmic microwave background 
(CMB) photons. As for heating, all species contribute to 
cooling, giving the total cooling rate 



For a single species, recombinations radiate the elec- 
tron energy ~ kT as photons at the rate Li = n e ni/3i(T). 
From Hn and Hem, the recom bination cooli ng coeffi- 
cients f3i(T) are, as on (T), from ISeatonl (|l959l). For the 
Hen radiative term we use the expression in iBlackl (|l98ll ) 
while we combine the approximation /3hoii = 3ryd qhoii 
i|Gould fc Thakurlll970l) with the second OHeii term listed in 
Table [1] (|Aldrovandi fc Peauignotlll973l ) for the dielectronic 
component. The total recombination cooling coefficients are 
provided in Table [T] 

For the cooling rate fro m collisi o nal ex citation of Hi we 
adopt the approximation of ISpitze r (1978): 



L eH = 7.3 x 10 



-32 T 3 

J m s 



-1 -118400/T 

n e nme 



(13) 



We note that for the temperatures relevant here, cooling 
losses due to collisional ionisation of hydrogen and collisional 
excitation and ionisation losses from helium are negligible. 

Finally, Co mpton scatter ing off CMB photons cools the 
gas at the rate (|Peebleslll968h 



Lc 



Aarak 



-T CMB (z) [T - T C mb(z 



(14) 



where <jt is the Thomson cross section, a is the radiation 
density constant, m e is the electron mass, c is the speed 
of light, and Tcmb is the temperature of the microwave 
background. 

In the current implementation, the simulations do not 
solve for the overlapping of ionisation fronts. Rather, they 
describe the passage of the first ionisation front across a neu- 
tral region. Reionisation will also be affected by the hydro- 
dynamical response of the gas, which is not included here. 
For low to moderately overdense systems, this has only a 
moderate effect on t he statistical properties of the resulting 
absorption systems (Me iksin fc Whitel l2001). The hydrody- 
namical effects are more important for reionisation in denser 
structures like minihaloes, which can be optically thick at 
the Lyman edge. Ionisation heating results in an overpres- 
sure that drives the gas out of the minihaloes, reducing their 
densities and making them less effective at slowing the ion- 
isation fronts. Estimates based on the photon consumption 
rate sug gest the role of the se systems in slowing the fronts 
is small jciardi et aTll2006t ). 

Another simplification is the absence of diffuse radi- 
ation. Radiative recombinations produce a diffuse radia- 
tion field throughout the ionised region. Besides the re- 
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combination rate, the intensity depends on the amount 
of clumping of the gas. Estimates range from a boost 
in the ionisation rate by an additional 10—40 per cent. 
jMeiksin fc Madaulll993l ; faaardt fc Madaulll996r i. After the 
gas has been ionised, the diffuse field may be accounted 
for by rescaling the overall radiation level, as the inten- 
sity of the radiation field has a negligible effect on the 
post-ionisation gas temperature. Radiative recombinations 
will also contribute to the reionisation process itself. The 
contribution, however, is generall y negligible for bo t h hy- 



drogen and helium reionisation (Shapiro & Giroux 



Meiksin fc Madau 



1987 



1993 



Miralda-Escude fc Ostrikeil Il99d ; 
Madau fc Meiksin! 11994). An exception is in dense regions 
in a scenario in which the Hem 1-front precedes the Hn I- 
front. In this case, assuming all photons produced by helium 
recombinations to the ground state of Hen are locally ab- 
sorbed by neutral hydrogen atoms, the time to ionise the Hi 
is 



"II 



5nl 



7.7 x 10 20 s 



(15) 



0.020 



1 + 5 



0- + z) 



where af e is the radiative recombination rate to the ground 
state of Hen, T4 is the gas temperature in units of 10 4 K, and 
8 is the local fractional overdensity. An allowance has also 
been made for 1.5 a dditional se condary electrons produced 
per photoionisation (|Shulllll979r i. Compared with a Hubble 
time of iHub ~ 14Gyr(l + z)~ 3//2 , the ionisation time goes 



tile 



^Hub 



! 640T4 



(l + z) : V2(i + ( 5)- 



(16) 



For the range of source turn-on redshifts considered (8 < 
z < 20), except in very overdense regions, this ratio will be 
large and hydrogen will remain largely neutral in the inter- 
val between helium and hydrogen ionisation. The effect on 
the temperature will therefore be small, except for very over- 
dense regions. After reionisation, the post-ionisation temper- 
ature will be accurately computed even in the highly over- 
dense regions since the temperature is set by the balance 
between ionisation heating and radiative c ooling with n o 
memory of the reionisation history retained l|Meiksinlll99l) . 

Heating by x-rays is automatically included, but not the 
effect of secondary electrons. The secondary electrons have 
a neg ligible effect on th e post-ionisation temperature of the 
gas i|Madau et alJll997h . Ahead of the I-front, however, the 
heating rates will be overestimated by a factor of 2-5, de- 
pending on optical depth to t he ionising photons (which de - 
termines their mean energies) (|Shull fc van Steenber Jl985h . 
While the tendency to approach thermal equilibrium reduces 
the impact of the secondaries on the final temperature, the 
gas temperature ahead of the I-front is somewhat reduced 
when the effect of the secondaries is included, as discussed 
in the next section. Because of the excessive computational 
expense involved in including the effect of secondaries, we 
neglect their role on the pre-ionised gas temperatures, but 
note the effect of secondaries must be included for applica- 
tions for which the pre-ionised gas temperatures are rele- 
vant, as f or instance the 21c m signature of the gas at high 
redshifts l|Madau et al.lll997f ). 

None of the above simplifications affects our findings 



presented in this paper, which concentrates on the post- 
ionisation thermodynamic properties of the IGM as it relates 
to the Lya forest. 



2.2 Implementation 

To solve the RT equations, we use the probabilistic formula- 
tion described in Abel et al. (1999), as extended to multiple 
species by Bolton et al. (2004) . In the probabilistic formula- 
tion, the ionisation rate per unit volume due to an incident 
photon flux per frequency interval, T „, is given by 



nF ■■ 



1 

dr 



T v a v ndrdv. 



is approximated as 



nV = 7T / ? > 



dv 



(17) 



(18) 



where t?, 



a v ndr is the optical depth of a single cell of 



width dr. This formulation conserves photons and does not 
lead to excessive ionisations through cells for which t^ c11 > 1. 
For multiple species, the probability of absorption is spread 
among the various specie s and the rat e of io nisation must 
be known for each species. iBolton et all (|2004 ) describes the 
split for absorption by Hi, Hei, and Hen. Using their short- 
hand p 1 — e~ T " and q x — 1 — e~ T " where the r"s are within 
the cell (as opposed to cumulative), the probabilities of ab- 
sorption by each species are 



Hi He I He II / -1 

= q p p I 1 — e 



He 1 Hi Hell 1 -1 
= q P P (1 



n HcII Hell Hi HcI / 

Pab S = q VP 1 1 



ID 
ID 
ID 



(19) 
(20) 
(21) 



„Hl„HcI„HcII 



„HcI„Hl„HcII 



,HoII„Hl„HoI 



-r^ 11 = t? 1 + t? c1 + t? cU . Equation [T7l with o v ndr replaced 
by the Pabs' s above is the form implemented in PMRT_mpi. 
We recast Eq. [TO] for each species in a similar manner, to 
prevent excessive heating. 

Both the PM and RT modules were coded in C and par- 
allelised using a message pas sing interface (MPI) (the PM 
component previously so fsee lMeiksin et al]|l999l )). MPI is 
suited to work on distributed memory systems, but works 
equally well on shared memory systems, albeit not with op- 
timal memory use. The RT module consumes most of the 
processing time. In the RT module, each LOS is processed 
by a single computer processor element (PE). In general 
the time spent by each process on a LOS scales linearly 
with the resolution along the LOS. In regions of high op- 
tical depth, refinements are generated which split the cell 
into low-optical-depth sections. The number of refinements 
is a factor in the PE load for a LOS. Load imbalance is 
possible because the time to process a LOS varies for each 
LOS, depending on the amount of structure. Load imbalance 
is alleviated by dynamic (first-come-first-serve) allocation of 
LOSs to processes. The master node performs the allocation 
and, consequently, has a significantly lower mean load. The 
load imbalance with the master node can be reduced by ei- 
ther 1) having one physical processor unit execute both the 
master and a slave process or 2) using many processor units, 
so that instead of one-of-few being underused, one-of-many 



Reionisation scenarios and the temperature of the IGM 5 



is underused and less time is spent in a non-load-balanced 
state. 

In addition to accurate modelling of the ionisation 
front, special attention is given to computing accurate post- 
photoionisation temperatures. This requires sufficient reso- 
lution in time and space to correctly model both the ion- 
isation structure of the ionisation front and the tempera- 
ture across it. This is particularly important when mod- 
elling reionisation by sources with hard photons like QSOs. 
The time step for the radiative transfer is computed, inde- 
pendent of the PM time step, from the cooling, ionisation, 
recombination, and density-change time scales. The mini- 
mum of these time scales is multiplied by a factor of 0.2, 
found through convergence tests using a uniform medium 
containing Hem and Hn fronts. Convergence is considered 
established if the maximum error in the temperature along 
the LOS (compared with a fiducial run using a factor of 0.1) 
is less than half a per cent. High spatial resolution at the 
front is guaranteed by the use of refinements. Any cell in a 
line of sight (LOS) which is optically thick at the Lyman 
edge (r > 1) is split into a sufficient number of slices such 
that r < 1 in each. Convergence tests were similarly used to 
establish this criterion by comparing with a fiducial run with 
r < 0.05 per slice. Refinements are normally only generated 
at the ionisation front. 

Integration over frequency of the ionisation (Eq. [2} 
and heating (Eq. I10p functions is performed using Gauss- 
Legendre quadrature over the intervals vm to unei and v>h b i 
to j^hcii and Gauss-Laguerre quadrature for the interval 
VHeii to oo where Vi is the ionisation threshold for species i. 
The integration parameters for each of these intervals was 
tailored to each of the input spectra using convergence tests 
similar to those described in the previous paragraph. 

Of course, the effort required to get the temperatures 
correct comes at the expense of computer cycles. More than 
20000 PE hours were required for the simulations presented 
here. PMRT_mpi was run in parallel on 8 or 16 PEs drawn 
from a dedicated cluster of IBM OpenPower 720s or a local 
cluster of Linux boxes of mixed type. 

Validation of the code was accomplished by the use of 
problems with known solutions and comparison with results 
produced by an independent code. The analytic solution to 
the position of a hydrogen ionisation front in a gas of uniform 
density irradiated by a monochromatic source with photon 
energ y hv = 13.6 eV (producing no heating) is easily deriv- 
able l lliev et "ai"1l2006l . for example). We simulated the pas- 



sage of an I-front through a gas with uh = 10 cm and 
T — 10 4 K ionised by a s ource producing p hotons at a rate 
of 5 x 10 48 s -1 (Test 1 of llliev et all (200g) ). Over a period 
of 500 Myr (a few times the recombination time scale), the 
error in the position of the front was never more than 6 per 
cent. There is no analytical solution for the gas tempera- 
ture if irradiated by a non-monochromatic source. However, 
other simulations exist to which we could compare. One of 
us (AM) has a radiative transfer code which implements 
an altogether different method with more physics than we 
included in PMRT_mpi (|Madau et al.lll997l '). For an a — 0.5 
power-law spectrum ionising a uniform IGM at the mean 
baryon density, this code (with the effect of secondary elec- 
trons switched off) produces post-ionisation gas tempera- 
tures of about 18,000 K to 20,000 K at z = 7. Including 
the effect of secondary electrons was found to reduce the 
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Figure 1. Source spectra. The ionisation thresholds for Hi, Hel, 
and Hell are indicated. 



post-ionisation temperature by about 200 K while the tem- 
perature in the gas ahead of the I-front, where helium ion- 
isation releases large numbers of electrons, was reduced by 
about 1000 K. Our PMRTjnpi code produces post-ionisation 
gas temperatures of 16,000 K to 23,000 K at this epoch. 
PMRT_mpi generates the correct ionisation structure and tem- 
peratures that agree with both expectations from observa- 
tions and the result of an independent code. We are confi- 
dent that given the physics included, the PMRTjnpi code is 
producing correct results. 



3 SIMULATIONS 

We use PMRT_mpi to simulate the passage of an ionisation 
front produced by sources with a variety of characteristic 
spectra. Each I-front computed may be considered to result 
from a single source with the adopted spectrum, or a group 
of several neighbouring sources with individual I-fronts that 
have already overlapped, producing a single advancing I- 
front through the IGM. Because we treat each line-of-sight 
separately, the resulting temperature distributions may also 
be considered as produced by distinct sources, and the re- 
sulting temperature statistics an ensemble average over the 
IGM. 

In the following, we describe the sources of ionising ra- 
diation ( i|3.ip . estimate the expected characteristic scales of 
their ionisation fronts at the reionisation epoch ( H3.2| l. de- 
scribe the volume of space simulated f i]3.3[) and demonstrate 
that the results have converged with resolution and are not 
dominated by cosmic variance ( H3.4p . 



3.1 Sources 

The source spectra (Fig. [1} were selected to emulate can- 
didate reionisation sources. For the power law, the lumi- 
nosity is given by L(u) rx where a — 0.5, which 
corresponds to a hard quasar/QSO/AGN sp ectrum. The 
miniquasar spectrum of iMadau et al.l (|200- 



L(x) 



+ 8oT 



1 < x < 20; P(x) 



is given by 
oc 8x _1 , x > 20 
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where x = v jvp. The starburst spectrum was produced by 
PEGAsf] i|Fioc fc Rocca-Volmerangelll997l 'l for a galaxy 30 
Myr after a burst of Pop III (zero metallicity) star forma- 
tion. The starburst spectrum has an effective spectral index 
of a e a = 7.4 just above the Lyman edge. The hybrid model 
begins with a starburst spectrum, and between z = 5 and 
z = 4 evolves into an a = 0.5 power law, mimicking a radi- 
ation field in which a QSO source dominates. 

The order of propagation of the Hn and Hem ionisa- 
tion fronts is different for the different sources. For a time- 
invariant spectrum, the Hem I- front precedes the Hn I-front 
if the spectrum is hard enough, specifically if a c ff < 1.8. 
The condition is met for both the power-law and miniquasar 
spectra and, indeed, the Hen was ionised prior to Hi in those 
simulations. The starburst spectrum certainly fails the con- 
dition as it has negligible intensity above vn<ni- The hybrid 
spectrum, in which a spectrum dominated by star bursts 
evolves into one dominated by a power law, leads to the Hn 
I-front preceding that of Hem. 

To ensure comparable results, all spectra were nor- 
malised to produce an incident hydrogen ionising flux of 
1.5 x 10 7 (1 + z) 2 s _1 m -2 . This corresponds to a typical 
flux level driving an ionisation front, as detailed in i|3.2l For 
example, it corresponds to a L[i/Hi] = 1-5 x 10 21 W Hz -1 
QSO source with a — 0.5 at a comoving distance of 5Mpc . 
The flux is adequate to ionise the full simulation volume by 
z — 3. We fix the incident flux level to a common level to 
ensure that differences in the post-reionisation temperature 
are due to the differences in the incident spectra rather than 
the time scale for reionisation. We also explore the effect of 
source turn-on redshift by performing, for each source, simu- 
lations with the source turning on at redshifts of z on = 8, 12, 
and 20. These were selected to cover the range of redshifts 
limited by the WMAP and QSO observations. Because of 
the common incident flux normalisation, these do not corre- 
spond to different reionisation scenarios in which the reion- 
isation of the Universe completes at different epochs, as the 
reionisation histories are nearly identical for z < 8. It is 
simply a device that allows us to ensure the results at late 
times are insensitive to the assumed turn-on redshift of the 
sources. 

The ionising radiation is projected as parallel rays nor- 
mal to the surface of the simulation volume. This configu- 
ration has the computational advantage of allowing all ra- 
diative calculations for a given ray to be done in a single 
column, independent of neighbouring columns. Since no in- 
formation about the thermal state of the gas is carried with 
the particles as they move from one cell to another, each col- 
umn can be interpreted as an individual line of sight which 
is independent of any other line of sight. We select 256 lines 
of sight arranged in a plane to assist in visually interpreting 
the passage of the front while providing enough lines of sight 
to deal with cosmic variance. The radiation field takes the 
plane-wave approximation and is not geometrically diluted 
by r -2 . A limited number of simulations performed with the 
radiation field geometrically diluted produced results similar 
to those without the dilution factor. 

For brevity, the labels listed in Table [3] will be used 



Table 3. Identifiers for the various models used in this paper. 



Model 


8 


^Oll 

12 


20 


Power Law 


PL08 


PL12 


PL20 


Miniquasar 


MQ08 


MQ12 


MQ20 


Starburst 


SB08 


SB12 


SB20 


Hybrid 


HY08 


HY12 


HY20 



henceforth to identify the various simulations with their 
combination of model spectrum and source turn-on redshift. 



3.2 Characteristic sizes of ionised regions 

A key effect for which radiative transfer is required is the 
modification of the source spectrum incident on a region by 
intervening gas, including shadowing by dense knots. The 
amount of modification depends on the amount and nature 
of the intervening gas. In this section, we estimate the typ- 
ical distances expected between a source and the I-front it 
produces, and show that our simulations match the distances 
for a wide range of plausible sources. This is important, be- 
cause the no-overlap approximation we make is only valid 
provided the I-fronts do not extend further than the char- 
acteristic size of an I-front at the time of complete reion- 
isation, when the porositj0 of the ionised gas in the IGM 
crosses unity. 

At the time of complete reionisation, the mean space 
density of the sources is the inverse of the characteristic 
volume of the ionised gas bubbles just before their com- 
plete percolation. For any given characteristic luminosity 
per source, the source space density may be estimated from 
the total emissivity of all the sources just prior to com- 
plete reionisation. We adopt an estimate for the emissiv- 
ity based on the constraints imposed by the measured IGM 
Lyq optical depth up t o z ~ 6 in a ACD M cosmology 
jMeiksin fc White! I2003L 12004 iMeiksinl I2005T ). At z = 6, 
the rate of Hi ionising photons per comovi ng volume wa s 
found to be A s Ri 4 x 10 50 a g 1 for 1 Mpc~ 3 (|Meiksinl 120051 ) , 
where as is the spectral index of the ionising sources at the 
Hi Lyman edge and h is defined by the Hubble parameter 
(Ho = 100/ikm s _1 Mpc _1 , where Ho is the present day 
value of the Hubble parameter) . For a source producing ion- 
ising photons at the rate iV 7 , the corresponding volume of 
ionised gas at the time of overlap is then Vi ~ N^/hs. The 
characteristic comoving size of an I-front at z — 6 is then 
(for h = 0.71), 



4-7T ns 



1/3 



9a 



1/3 



N-. 



10 54 ! 



1/3 



Mpc. 



(22) 



This scale agrees with estimates made by others. 
IWvithe fc Loebf(|2004h determined at overlap the ionised re- 
gions are constrained to have scales of 5 to 50 ft -1 Mpc 



1 http://www2.iap.fr/users/fioc/PEGASE.html 



2 The porosity of the ionised gas is the product of the source 
number density and ionisation volume produced by an individual 



source. 
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comoving. Similar scales have bee n inferred using semi- 
analytical reionisation approaches dFurlanetto fc Ohl [20051 ; 
iFurlanetto et al.ll2006l : ICohn fc Chandl2007i ). 

In the following sections, we estimate the characteris- 
tic I-front size at the time of percolation for Lyman-break 
galaxies, low luminosity AGN, and a few other plausible 
sources. 



3.2.1 Lyman-break galaxies 

Lyman-break galaxies discovered at z > 6 have been sug- 
gested as principal sources of the reionising photons, de- 
pending on the uncertain escape fraction of ionising radi- 
ation from the galaxies, the nature of the dominant pop- 
ulation of stars in the galaxies, and the clumping fac- 
tor of the IGM dBunker et alj 12004 IStiavelli et all |2004 
lYan fc Windhorst! l2004r i~ The sources have characteristic 
UV (1500 A) luminosities of L v « 10 28 - 10 29 ergs" 1 Hz" 1 
ISteidel et al.lll999l ; [Bunker et al.ll2006l ). The corresponding 
rate of ionising photon production depends on the conver- 
sion rate between ionising photons and UV light. Using the 
values for Population III, Po pulation II and Pop ulation I 
(solar abundance) stars (e.g. IStiavelli et al.l 12004 ). we ob- 
tain Nj = (4.8, 1.0 , 0.17) x 10 55 s _1 , respectively. Allow- 
ing for an escape fraction / csc , Eq. [22] gives for the char- 
acteristic comoving I-front radius at the t ime of overla p 
Ri ~ 15(/ csc /0.1) 1/3 Mpc for a s ^2- 7.4 (|Meiksinll2005l) . 
The uncertainty in the ionising photon conversion rate and 
escape fraction corresponds to an additional factor of un- 
certainty of about two in the size. Addressing the size of 
the ionised bubbles in which QSOs reside prior to turn-on, 
lAlvarez fc Abell (|2007l ) show that galaxies ionised regions on 
the order of tens of comoving Mpc shortly before overlap. 

It is possible the IGM was reionised by even rarer, more 
massive and more luminous systems, su ch as massive post- 
starburst galaxies l|Panagia et al.1 12005). Allowing for Pop- 
ulation III stars and a high escape fraction, typical I-fronts 
at the time of overlap could approach comoving radii of 
~ 30 - 50 Mpc. 

3.2.2 Active Galactic Nuclei 

The number counts of luminous QSOs fall short of the 
amou n t required to reionise the IGM (|Yan fc Windhorst] 
12004 iMeiksinl 120051 ). Low luminosity AGN, however, 
are plausible sources of non-stellar ionising radiation 
(Rico ttTfc Qstrikerl20 04f). A population of galaxies harbour- 
ing 10 6,5 Mq black holes shining at their Eddington luminosi- 
ties would require a space d ensity of n a ~ 2 x 10 _4 Mpc -3 to 
reionise the IGM by z = 6 fMciksin 2005|). The correspond- 
ing characteristic comoving I-front radius at the reionisation 
epoch is then Ri » n g ~ 17 Mpc, comparable to that ex- 
pected for Lyman-break galaxies. 

It is likely that the sources that ionised Hen to Hem 
were QSOs, as few other objects produce an adequate sup- 
ply of sufficiently hard photons. If the QSOs were low lu- 
minosity QSOs, as above, then the above estimate for the 
Hn fronts would apply to the Hem ionising fronts as well. 
If the rarer luminous QSOs dominated the reionisation of 
Hen, the characteristic I-fronts at the time of overlap would 
be even larger. For a comoving space density at 3 < z < 5 



of 10~ 7 - 10" 6 Mpc" 3 l|Meiksinll20"05l ). the characteristic co- 
moving sizes would be 100 — 200 Mpc. 

3.2.3 Other sources 

It is possible that more common structures ionised the IGM, 
such as collapsed systems smaller than Lyman-break galax- 
ies which mer ged into larger systems. In the simulations of 
iGnedinl ([2000), reionisation is dominated by systems with 
total stellar masses exceeding 3 x 10 s Mq . As only a few such 
systems are responsible for ionising the gas in a 4/i _1 Mpc 
comoving volume, the characteristic I-front size at overlap is 
a few comoving Mpc. Our models are just marginally appli- 
cable to such sources, but more relevant to the rarer larger 
I-fronts that would occur prior to overlap than the smaller 
ionised regions. 

In the miniquasar model of lKuhlen fc Madaul (|2005l ). an 
individual source does not ionise hydrogen much beyond a 
comoving distance of 20 kpc. In this case, the characteristic 
I-fronts are much smaller than those we model. Our mini- 
quasar source spectrum correspo nds to luminositie s muc h 
larger than those considered by IKuhlen fc Madaul (2005). 
Nevertheless, we adopt the miniquasar model as a physi- 
cally motivated alternative hard spectral shape. As we show 
later, the results are nearly indistinguishable from the pure 
power-law case. 

A more exotic source of ionising radiation is some form 
of decaying particle. Various hypotheses have been sug- 
gested. Our results are not relevant to a scenario in which 
such particles were responsible for the reionisation of the 
IGM, as the sources are completely localised and ubiqui- 
tous. 

3.3 Simulation volume 

All simulations were performed in a (25/i _1 Mpc) 3 comov- 
ing volume. This is large enough to include the typical pre- 
overlap I-front for t he models discussed in ^3.21 A ACDM 
model was assumed l|Spergel et al.|[2003l ). with parameters: 
h = 0.71, n b h 2 = 0.022, fi m = 0.268, fi„ = 0.732, where 
fib, fi m , and fi„ have their usual meanings as the contribu- 
tions to fi from the gas, all matter, and the vacuum energy, 
respectively. For the hydrogen fraction by mass, Y = 0.235. 

The initial density perturbations were created by dis- 
placing a uniform grid using the Zel'dovich approximation. 
The initial power spectrum of the density fields was a COBE- 
normalised power law with index n — 0.97. The same initial 
conditions were used for all simulations, except for the con- 
vergence and cosmic variance tests. Since there is no feed- 
back from the RT to the PM code, all the runs have identical 
gas densities. The simulations were evolved to a redshift of 
3. 

3.4 Cosmic variance, convergence 

The results are subject to cosmic variance and bias from the 
resolution of the A-body code and the RT grid, but not sub- 
stantially. In this section we describe two extra simulations 
which ascertain the effects of cosmic variance and a change 
in the resolution parameters by comparing these with the 
PL08 run which they most closely match. We use the p — T 
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Figure 2. Convergence with resolution of the TV-body code and 
the effect of cosmic variance. The top panel illustrates the p — T 
distribution for the ionised gas for two N p = 256 3 , N g = 256 reali- 
sations with the same power law model (black/solid and red/dash- 
dot) and a N p = 512 3 , N g = 512 realisation (purple/dashed). 
The contours are at dP/dlog (p/(p>) /dT = 3.6 X 10 _5 K~ 1 and 
3.6 X 10 — 4 K — 1 . The bottom panel shows the mass- weighted tem- 
perature distributions for the ionised gas for the same models. 



and temperature distributions (Fig. [2] top and bottom pan- 
els) to illustrate the effects. The extra simulations also show 
that the lack of advection in the simulations is not a source 
of significant error. 

There are two parameters that control the spatial reso- 
lution of the simulation: the number of particles in the PM 
simulation, N p , and the mesh size of the gas density grid, 
N g . To estimate the variation due simply to cosmic variance, 
a simulation was run with a different realisation of the initial 
density fluctuations at the same resolution as the main body 
of runs (N p = 256 3 , JV 9 = 256). In Fig.H the difference be- 
tween the solid (black) and dot-dashed lines (red) indicates 
the cosmic variance. A simulation with 8 times the mass 
resolution (N p = 512 3 ) and twice the RT mesh resolution 
(N g — 512) produced a qualitatively equivalent distribution 
(Fig. [2] dashed line), bracketed by the two lower resolution 
realisations. 

The current code does not advect thermodynamic quan- 
tities. Hence if a dense halo has a high velocity, the gas in 
the halo is not properly modelled. The error from the lack 
of advection is not easy to quantify, but it is not believed to 



be large. The error is more prominent when cells are smaller 
since moving gas is more likely to cross a smaller cell. The 
comparison of the 256 3 runs to the higher-resolution 512 3 
run (Fig. [2} demonstrates little variation. 



4 SIMULATION RESULTS 

The simulations provide information for two epochs of inter- 
est: the period during which the gas is being ionised and the 
post-ionisation epoch. Our primary focus in this paper is the 
post-ionisation temperature of the gas, but future telescopes 
such as the Low Frequency Array □ (LOFAR), the Mileura 
Widefield Array^ (MWA), and eventually the Square Kilo- 
metre Array^ (SKA), will be able to probe the epoch of ioni- 
sati on itself via the red shift ed 21cm line of neutral hydrogen 
(see iFurlanetto et alj [20061 . for a review). The reionisation 
epoch will be the topic of a future paper. 

In this section, we focus on the thermal state of the 
gas as established by the passage of an individual ionisation 
front. The reionisation process involves the overlap of such 
fronts. The ionisation fractions of the gas after the reionisa- 
tion epoch will be constantly reset as the gas sees ionising 
photons from an increasing number of sources, both because 
of the time needed for photons to reach the gas and as new 
sources turn on (while older ones die). These effects are not 
accounted for by our single I-front models. The temperature 
of the gas, however, is nearly insensitive to the continual 
resetting of the level and shape of the ionisin g photon back- 
ground ijMeiksinl 1 19941 : iHui fc Gnedlnl 1 19971 ). It is instead 
determined by the reionisation process itself and the subse- 
quent evolution of the physical state of the IGM. We here 
explore the evolution of IGM temperature as computed in 
our simulations. A principal goal is to determine if the tem- 
perature of the ionised IGM is dependent on the nature of 
the ionising source. We find that it is, and in £14.31 we quan- 
tify and discuss the magnitude of this effect and some of its 
consequences for the ionisation structure of the IGM. 

The simulation results also permit us to explore a few 
other issues related to the reionisation of the IGM. Clumping 
of the gas will generally impede the propagation of ionisa- 
tion fronts. We evaluate the importance of this effect in i )4.2l 
Because of the inclusion of helium ionisation in our simula- 
tions, we are able to examine the ratio of Hen and Hem to 
Hi in the ionised IGM prior to the epoch of complete Hen 
reionisation. We discuss these results in ij4.5l 

But first, to justify the effort put into the implemen- 
tation of RT in our code, we make a comparison with the 
optically thin approximation in £|4. 1 1 



4.1 



Comparison with the optically thin 
approximat ion 



Neglecting RT can underestimate temperatures by a fac- 
tor of a few (lAbel fc Haehnehj 1 19991 : IBolton et all 120041 : 
iBolton fc HaehneliJ 12007^ 7 The importance of including ra- 
diative transfer effects during reionisation to obtain accu- 
rate temperatures is best illustrated by comparing against 

3 www.lofar.org 

4 web. haystack. mit.edu/ arrays/M WA / site / index. html 

5 www.skatelescope.org 
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Figure 3. Effects of optically thin approximation. The top panels 
illustrate the p— T distribution without (left) and using (right) the 
optically thin approximation. The lower panels show the modifica- 
tion to the volume- weighted (left) and mass-weighted (right) tem- 
perature distributions produced without (solid) and with (dot- 
dashed) the approximation. The volume-weighted distribution is 
affected more because most of the volume of the Universe contains 
underdense regions, which retain a memory of the reionisation 
process. 
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Figure 5. Fraction of the simulation volume reionised for three 
models, compared with the case for a uniform density IGM at the 
mean density. The hybrid model follows the same history as the 
starburst model prior to z = 6. A simple model in which a uni- 
form density gas (p = Q g p c [l + z] 3 ) is ionised by a constant flux 
of ionising photons is plotted as the solid line. There are three 
separate redshifts at which the ionising flux is turned on. All the 
models initially ionise the gas faster than the uniform density case 
as the I-front sweeps into underdense regions. Clumping eventu- 
ally slows down the growth of the Hn filling factor by an amount 
depending on the spectral shape. 

density is produced in both simulations, resulting from the 
establishment of thermal b alance between photoionisation 
heating and atomic cooling |Meiksinlll99i ). 
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a simulation using the optically thin approximation (OTA). 
Simply put, the OTA means all points not near a source 
receive the same average radiation field. The approximation 
works fairly well because any gas dense enough to self-shield 
is also dense enough to reach thermal balance. So the post- 
ionisation temperature is essentially insensitive to the details 
of reionisation. The approximation fails, however, to prop- 
erly treat low density gas, because the time to reach thermal 
equilibrium in low density gas exceeds the Hubble time at 
high redshifts. As a consequence, the low density gas retains 
a memory of the reionisation details. This is particularly im- 
portant when low density gas is shielded from the oncoming 
I-front by dense structures because of their role in hardening 
the spectrum of the radiation field. 

To compare our results with those using the OTA, 
we performed a simulation in which the cumulative opti- 
cal depth to any point was set to zero, mimicking the OTA. 
The p — T and temperature distributions for the PL08 model 
with and without the approximation are illustrated in Fig. [3] 
Overall the distributions are similar. Since different regions 
are exposed to different spectra without the OTA, there is 
more spread than when the OTA is used. The hardening of 
the spectrum due to selective absorption of low-energy pho- 
tons heats the gas, particularly in the less dense regions. A 
high density spur of decreasing temperature with increasing 



4.2 Clumping 

It has long been recognized that clumping of the 
IGM gas may substantially slow the propagation of I- 
fronts d ue to the increased ra t e of radiative recombi- 
nations dShapiro fc Girouxl Il987l ; iMeiksin fc Madaul Il993l ; 
iMadau fc Meiksinl Il994h . Estimates of the importance of 
clumping have varied, but rece ntly tend toward only a mod- 
erate slowdown of the I- fronts (Sokasian ct aL 2003; Meiksinl 
I2005I ; ICiardi et alj|2006h wit h the resistance in creasing with 
time as clumping increases (jlliev et al.ll2005l ). In our sim- 
ulations, the slowing of the I-fronts in individual lines of 
sight leads to shadows in the ionisation maps, as shown in 
Figure |4] 

We have estimated the role clumping may play in delay- 
ing reionisation by comparing the growth of the H11 filling 
factor in the simulations with that predicted for a uniform 
IGM at the mean baryon density up to z — 6, at which point 
overlapping I-fronts would complete the reionisation. Intro- 
ducing structure into the IGM actually results in a moderate 
increase in the rate at which the filling factor of ionised gas 
grows in the early stages, as shown in Fig. [5] This is because 
most of the volume of the Universe is underdense. Once half 
of the volume is reionised, the filling factor converges to the 
uniform density prediction, with a small slowdown depend- 
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Figure 4. The Hi fraction at z = 6. From top-left to bottom-right the panels correspond to PL20, MQ20, SB20, and HY20. (The hard 
source has not yet turned on in the hybrid model, so the simulation at this point is identical to the starburst run.) The ionising flux is 
incident from the left. Distances are in proper units. 



ing on the spectral shape of the source. The evolution of the 
mass- weighted fraction of ionised hydrogen follows a similar 
trend including the more rapid rise than the uniform den- 
sity case at early times. By z — 6, however, the fraction 
grows substantially less rapidly, with 50 per cent to a fac- 
tor of two less mass ionised than for the uniform density 
prediction. The slower growth for the mass-weighted case vs 
the volume filling factor is expected since it takes longer for 
the I-front to penetrate the densest regions which contain 
proportionally more mass. 

It is possible the amount by which the I-front slows 
down is underestimated because of the deficit of small dense 
structures like Lyman Limit Systems in iV-bod y simulations 
l|Gardner et al.lll997l ; iMeiksin fc White] I2004T I. A definitive 
result may need to await hydrodynamical simulations that 



reproduce the statistics of Lyman Limit Systems and denser 
intergalactic structures. 



4.3 Temperature 

Figure [6] is a map of the temperature distribution at a red- 
shift of z = 6. Figure [7] shows a similar map for the power- 
law model at z = 3. In both maps, only regions in which 
the Hi was ionised (/hi < 0.1) at z = 6 are shown. (The re- 
maining gas would have been ionised by a different ionisation 
front or fronts by this time.) We immediately see a variety 
of effects that will be explored quantitatively later. First, 
the gas temperature in the miniquasar model is virtually 
identical to that of the simple a = 0.5 power law. Second, 
the starburst model produces the coolest ionised gas. (At 
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Figure 6. The temperature at z = 6. From top-left to bottom-right the panels correspond to PL20, MQ20, SB20, and HY20. (The hard 
source has not yet turned on in the hybrid model, so the simulation at this point is identical to the starburst run.) The ionising flux 
is incident from the left. Only regions for which /hi < 0.1 are shown. The remaining gas (masked as blue) would be ionised by other 



this stage, the hybrid and pure starburst models are iden- 
tical.) Third, the highest temperatures are typically found 
just behind the Hn I- front. Fourth, a "streaking" effect is ap- 
parent. The gas temperature downstream of a dense clump 
of gas is enhanced compared with the surrounding g 
result of delayed ionisation. The enhancement persists un- 
til z — 3. Finally, when compared with the gas density, all 
the models produce gas temperature structure that traces 
the large-scale gas structure, as illustrated in Fig. [7] The 
relation, however, is not one of simple proportionality, as 
we will see below. For instance, because the most recently 
ionised gas tends to be hotter at a given density, the tem- 
perature tends to increase towards the right (because the 
I- front passes from left to right), as shown in Fig. [6] 

The gas temperatures distinguish the various models 



at all redshifts. The temperature distributions of the gas 
ionised by z = 6, both volume-weighted (Fig. [8} and mass- 
weighted (Fig. [9]), show clear differences for the different 
model spectra. At z — 3, for the power-law and miniquasar 
models, the temperatures span 8 to 28 x 10 3 K (90 per cent of 
the gas mass with 5 per cent below the range and 5 per cent 
above) with a mass- weighted mean of 15 x 10 3 K . The hybrid 
model has a hotter tail in its distribution, ranging from ~ 10 
to 31 x 10 3 K with a mass-weighted mean of 17 x 10 3 K . 
Gas two to three times cooler is produced by the starburst 
model. For this model, the temperatures at z — 3 range from 
2000 to 17000 K with the mass- weighted mean at 9000 K . 
The distributions are, however, highly skewed with mass- 
weighted modes at 11000 K for the power-law, miniquasar, 
and hybrid spectra and 3000 K for the starburst spectrum. 
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Figure 7. Density distribution (left) and temperature map at z = 3. The temperature map is from PL20. Only regions for which the 
gas is ionised prior to z = 6 are shown, the remaining regions are masked as blue. 
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Figure 8. The volume-weighted temperature distribution at z = 
6, 5, 4 and 3 for the gas ionised before z = 6 only. From top-left to 
bottom-right the panels correspond to the power-law, miniquasar, 
starburst, and hybrid model. The colours/line type distinguish 
the rcdshifts, with 2 = 6 black/dash-dot-dot-dot, 2 = 5 red/dash- 
dot, 2 = 4 blue/dashed, and 2 = 3 green/solid. Note the change 
in the temperature scale for the various models. 
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Figure 9. The mass-weighted temperature distribution at z = 6, 
5, 4 and 3 for the gas ionised before 2 = 6 only. From top-left to 
bottom-right the panels correspond to the power-law, miniquasar, 
starburst, and hybrid model. The colours/line type distinguish 
the redshift, with 2 = 6 black/dash-dot-dot-dot, 2 = 5 red/dash- 
dot, 2 = 4 blue/dashed, and 2 = 3 green/solid. Note the change 
in the temperature scale for the various models. 



There are few direct determinations of the IGM gas 
temperature. The measured Doppler parameters on their 
own provide only an upper limit since the lines may be veloc- 
ity broadened (e.g. due to microturbulence) . If metal lines 
are present, however, the thermal and kinematic contribu- 
tions are s eparable. By s imult aneously fitting Civ and Siiv 
absorbers, iRauch et al.l |l996) found average temperatures 



of ~ 38000 K for systems with log 10 nm > 14. Only the hy- 
brid model is able to achieve temperatures of T ~ 40000 K 
(Figures [8] and O. 

We have tested that the redshift at which the source 
turns on is not a significant factor. Save for a slight shift to 
higher temperatures for the z on — 8 models, the curves are 
nearly identical. This is because the incident ionising flux 
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Figure 10. The p—T distribution at z = 3 for the gas ionised by 
2 = 6 only. From top-left to bottom-right the panels correspond to 
PL20, MQ20, SB20, and HY20. The contour intervals are spaced 
by a factor of vTO. 



has been normalised to a common value for all the cases so 
that by z £ 8 the reionisation proceeds in a nearly identical 
manner. Changing the epoch of reionisation would change 
the temperature in the low density gas at later times. This 
effect is not explored here. 

As we showed above, the use of full radiative transfer 
instead of an optically thin approximation not only further 
heats the gas, but it spreads the p — T distribution away 
from a tight power law (Fig. [3}. The spread in temperatures 
is greatest for gas with p/p < 1. The spectrum of the ion- 
ising radiation is modified by preferential absorption of the 
lower frequency photons as they pass through the gas. The 
modification hardens any spectrum, but is most influential 
to spectra that already have a hard component. Hence, the 
spread should be largest for harder spectra. Figure [TOl which 
maps the p — T distribution at z = 3 for the ionised gas for 
PL20, MQ20, SB20, and HY20, confirms the larger spread 
in the harder spectra. Also confirmed is the similarity of 
the state of the gas in the models with the power-law and 
miniquasar spectra. 

Pseudo-hydrodynamical models of the IGM usually 
adopt a polytropic equation of state for the gas. We fit a 
polytropic relation, T = To(p/p) 7_1 , to the p — T distribu- 
tions at z = 6, 5, 4 and 3. The results are shown in Ta- 
ble |4] Note that the errors are for the coefficients and not 
an indication of the spread of the p — T distribution, which 
is not particularly well-described as a single polytropic de- 
pendence. For comparison, for the OTA run, for which local 
heating is balanced by local cooling, we find To = 13780 ±40 
and 7 = 1.527 ± 0.003 at z = 3, with less spread than found 
when radiative transfer is incorporated. Together, these very 
different simulations imply a local balance between heating 
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Figure 11. Distribution functions for the ionisation rates at 
2 = 6. The distributions from the various species are distin- 
guished as indicated. Only rates from regions with ionised hy- 
drogen (/hoII < 0.1) are shown. From top to bottom, the panels 
correspond to the power-law, miniquasar, and starburst models. 
(The hybrid model is identical to the starburst model at this 
epoch.) The only contribution to the starburst distribution comes 
from hydrogen. Note the different T range for the starburst model. 
Also note that the distributions for the individual species are nor- 
malised independently for clarity. All the distributions show ex- 
tensive tails toward low values resulting from the absorption of 
photoionising photons by the IGM. 



and cooling still dominates the setting of the polytropic in- 
dex. 

Generally 7 = 1 is assumed to be the theoretical lower 
limit, corresponding to isothermality. Most of the fits indeed 
show 7 > 1. The exception is the value of 7 for the hy- 
brid model. Soon after the QSO starts ionising Hen, lower 
values of 7 occur, including 7 < at z — 4. This arises 
because of the poor representation a polytrope gives of the 
wide spread in temperatures. Any applications restricting 
to 7 > 1 may therefore be excluding a description of the 
thermal behaviour of the IGM during the Hen reionisation 
epoch. At the other extreme, for all models the index is lower 
than the adiabatic index, 7 a( j = 5/3, reflecting the influence 
of radiative cooling in high-density regions. 

The polytropic fits resul ting from the s imula tions may 
be compared with those of [Schave et al.l (|2000l ) to Keck 
HIRES spectra. These authors find T « 2.2 ± 0.2 x 10 4 K 
and 7«1.0±0.1atz~3, values most consistent with the 
late Hen reionisation hybrid model. 
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Table 4. Fits to the polytropic relation, T = Tb(p/p) 7 1 of the gas for which /hi < 0.1 by z = 6. The errors are to the fits and not 
indicative of the distribution in the p — T plane. 



z 







6 


5 


4 


3 


PL20 


Tb[K] 


17900 ± 100 


16900 ± 150 


15400 ± 140 


14600 ± 200 


7 


1.29 ± 0.01 


1.40 ±0.01 


1.52 ±0.01 


1.58 ± 0.01 


MQ20 


T [K] 


17900 ± 120 


16950 ± 170 


15500 ± 150 


14500 ± 160 


7 


1.34 ±0.01 


1.44 ±0.01 


1.52 ±0.01 


1.57 ± 0.01 


SB20 


2b [K] 


4200 ± 130 


4070 ± 50 


3720 ± 20 


3515 ± 11 


7 


1.44 ±0.03 


1.43 ±0.01 


1.47 ± 0.01 


1.553 ± 0.002 


HY20 


2b IK] 


4200 ± 110 


4040 ± 50 


29250 ± 250 


18600 ± 100 


7 


1.42 ± 0.03 


1.42 ±0.01 


-0.0363 ± 0.009 


1.134 ± 0.005 



4.4 Ionisation rates 

The varying optical depth to ionising photons results in fluc- 
tuations in the ionising background. The ionisation back- 
ground is parametrized by the ionisation rate, T, which is 
the number of photoionisations per atom per unit time. A 
long tail towards low ionisation rates is found, as shown in 
Fig. [11] The sharp cut-off at the high end corresponds to low 
optical depth to the source leading to negligible filtering of 
the source spectrum. Filtering decreases the ionising photon 
flux, lowering T. 

A tail in the r-distributiqn is sim ilarly found in the sim- 
ulations of Masclli & Ferraral l|2005l) . although without the 
sharp cut-off at the high end. They model a uniform source 
distribution filtered through the IGM, accounting for the ra- 
diative transfer of the ionising photons using a Monte Carlo 
algorithm. Th e absence of the s h arp u pper cut-off in the 
simulations of iMaselli fc Ferraral (2005) is consistent with 
the discreteness effects of their Monte Carlo radiative trans- 
fer approach. In the presence of a population of randomly 
distributed local ionisation sources, iMeiksin fc White! (2003) 
show that the sharp cut-off would be replaced by a power- 
law tail varying as T^j?' 5 . 



4.5 Helium 

The presence of helium alters the temperature of the IGM 
after reionisation by an amount that depends on the reion- 
isation scenario. The epoch of helium reionisation (ionisa- 
tion of Hen to He in) is still unknown. Measurements of 
the Hen Lya optical depth suggests i t occurred at 



l|Zheng et a l. 2004; Rei mers et al. 12005). which is consistent 
with the expected epoch of Hen reionisation by QSO sources 
with soft spectra fMeiksin 2005 1 ) ■ Measurements of the Ly af- 
forest Dopple r parameter, 6 , perm it an estimate of the IGM 
temperature; iTheuns et~aT] i|2002l ') and lRicotti et all (|2000l ) 
claim a temperature jump of about a factor of two to three 
at z > 3. During the reionisation process, and prior to its 
completion, large fluctuations in the Hen to Hi absorption 
signatures may be expected, as the Hen- ionising metagalac- 
tic UV background will show large spatial variations. 

We examine the predictions for these fluctuations from 
our simulations prior to the completion of Hen reionisation, 
by which time the Hem I- fronts have completely overlapped. 
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Figure 12. The importance of helium in the hybrid model. Illus- 
trated is the evolution of the temperature of a single line of sight 
in a HY08 simulation with (top) and without (bottom) helium. 
The transition between starburst and power-law spectrum occurs 
between z = 5 and 4. 



4-5.1 Thermal effects of helium 

The highest temperatures are produced by the hybrid model. 
The high temperatures are directly attributable to the pres- 
ence of helium. In the hybrid model, the irradiating spec- 
trum undergoes a transition from a starburst to a power-law 
spectrum between z — 5 and 4. As a consequence, the hy- 
drogen I-front precedes the helium I-front. Since the temper- 
ature of ionised gas cannot be raised significantly by chang- 
ing the radiation field, no matter how hard the spectrum 
(barring extreme cases like a pure x-ray spectrum), without 
helium the harder power-law spectrum has no effect on the 
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temperature of the gas previously ionised by the starburst 
spectrum. Gas that has not yet been ionised is heated to 
higher temperatures than that ionised by the starburst spec- 
trum since more energy, h(u — vo), is liberated per ionisation. 
In Fig. 1121 bottom panel, the differential heating in a hybrid 
model without helium (Y — 0) is illustrated. The gas ionised 
prior to the transition to the power-law spectrum remains 
unaffected by the transition. The gas previously unionised 
is ionised and heated to higher temperatures. With helium 
(top panel of Fig. I12[) the transition of the spectrum heats 
all the gas along the line of sight. This is not surprising since 
the ionisation of helium introduces additional heating. What 
is notable is that the gas is heated to higher temperatures 
than for the case of a power-law spectrum ionising a com- 
pletely neutral medium (~ 30000 K versus ~ 16000 K , see 
Fig. I) . 

A He in I-front passing through Hn-dominated gas leads 
to a larger jump in temperature than a Hn front passing 
through He in-dominated gas. (Note that in both cases, the 
final state is fully ionised.) The temperature-entropy rela- 
tion, Eq. [5] accounts for the differential temperature jumps. 
In the case of the Hen-Hem transition in a Hn-dominated 
gas, the mean molecular weight /x is reduced by only 7 per 
cent, leaving any gain in entropy to translate into a gain in 
temperature. In the alternate case of the Hi-Hii transition 
in a gas where Hem is the dominate helium species, /j, drops 
by 45 per cent, almost halving the gain in entropy. There is 
still a net gain in temperature, but it is only about 10 per 
cent. 

The shallow p — T profile in the hybrid model results 
from hotter low-density gas compared with the other models 
(Fig. 1101 lower-right panel). The high temperatures following 
the recent Hen reionisation are retained in the low-density 
gas due to the long time for thermal equil ibrium to be es- 
tablished in underdense gas (|Meiksinlll994l ). 

4-5.2 Correction to the Hi abundances at z < 6 

Prior to overlap of the Hn I- fronts at z ~ 6, the ionisation 
state of hydrogen is properly modelled by a single source, as 
we have in our simulation. Equivalently, a single source is 
sufficient for modelling the helium fractions prior to overlap 
of the Hen I-fronts at z ~ 3. However, at z < 6, a single 
source is insufficient to fix the hydrogen fractions, since once 
the Universe is reionised a given region becomes exposed to a 
large number of ionising sources. In particular, low-density 
regions once shadowed by Lyman-limit systems along the 
line of sight to the source driving an advancing I-front would 
be swept over by other I-fronts after the reionisation epoch. 

In order to compare the helium ionisation fractions with 
the hydrogen, we need to remove the effect of the Lyman- 
edge shadows on the Hi from the simulation data. We do so 
by correcting the Hi abundances by setting Tui to a uniform 
value at any given redshift and solving Eq. [TJfor nm in the 
static case. Estimates of the expected distribution of Thi 
after the reionisation epoch suggest it i s nar rowly peaked, 
especially for z < 4 (|Meiksin fe White! l2003t ). We set Thi 
to the mean value for the unshadowed gas at the given red- 
shift[f| A negligible change to the temperatures will be gen- 

6 This choice is arbitrary. We could have assigned the expected 



erated by an increased ionisation rate after other sources are 
revealed, so that the local temperature produced during the 
reionisation phase will still apply and so sets the value of 
the radiative recombination coefficient ami in Eq. [1] 

The solution to nm is still not entirely correct, since 
the dense Lyman-limit systems should contain self-shielded 
regions within which the correction eliminates. The filling 
factor of the missing self-shielded regions in the corrected 
version, however, is much less than that of the shadowed 
low-density regions. We therefore applied the correction to 
all the Hi data at z < 6. 



4-5.3 Comparison of Hen vs Hi absorption prior to 
complete helium reionisation 

By having two ionisation states, both at higher energies than 
that of Hi, helium provides a means of obtaining informa- 
tion about the abundance of high-energy photons. The ratio 
of the hydrogen and helium column densities is a function 
of the spectral shape of the ionising radiation. Coinciding 
Hi and Hen absorption features in QSO spectra ha v e bee n 
compared bv lZheng et al.l i|2004l ) and lReimers et all l|2005h . 
who find large fluctuations in the ratio of Hen to Hi column 
densities at z £ 3. The presence of fluc tuations plac e s con - 
straints on helium reionisation models. Ideser et all (|2005l ) 
find the patchiness requires short QSO lifetimes (< lOMyr) 
in models which attribute the pat chiness to the discre teness 
of the QSO spatial distribution. iBolton et al.l (|2006h have 
argued for a model in which the relative number of rare 
Hen-ionising QSOs compared with abundant He i- ionising 
star-forming galaxies sets the fluctuation distribution. 

We confine our estimates to the amount of fluctuations 
expected behind Hem I-fronts prior to their complete over- 
lap, i.e., before the epoch of helium reionisation has com- 
pleted. We base the analysis on the simulation results at 
z — 3 (Fig. 1 13 p . which is when observations suggest Hen 
reionisation was nearing completion. Local sources of Hi- 
ionising photons could introduce further fluctuations than 
those we find. Our simulations thus only set a lower limit to 
the level of fluctuations expected, arising principally from 
shadowing and attendant fluctuations in the Hen ionisation 
rate, during the final stages of Hen reionisation. 

The hardness index r\ = nHeii/nm is a sensitive probe 
to the shape of the ionising background. Ionisation fraction 
fluctuations affect the spread in n values while the mean 
value is set by the mean spectral hardness. Here we con- 
centrate on the relative spread of the fluctuations in n, as 
this is independent of the shape of the spectrum. We do re- 
fer to definite values for r\ from our simulations rather than 
arbitrarily rescaling them, but the physical effects we shall 
describe are not specific to any specific value n. 

The ionisation states of hydrogen and helium are ex- 
pected to be linearly correlated, but dependent on the hard- 
ness of the local radiation field. This follows from radiative 
equilibrium, for which the rate equations for hydrogen and 



value at the corresponding redshift, as determ ined by matching 
to th e measured Lya effective optical depth iMeiksin &; White] 
2004), but since n is proportional to the ratio THi/rHeiii it is 
anyway fixed only up to an overall rescaling factor. 
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Figure 13. The Hen fraction at z = 3. From top-left to bottom-right the panels correspond to PL20, MQ20, SB20, and HY20. The 
regions in which hydrogen is not ionised by z = 6 are shaded blue. 



helium are derivable by setting hi = in Eq. [T] 

ninTHi = n e nHiia?Hii (23) 

^HoIirHoII = W e nHeIIiaHeIII- 

For the highly ionised gas in the simulations, /hiThi — 
n e QHn and /HenTHeii — n e a H ciii, giving 

/Hell QHelll Th I 



/HI 



oh ii Thcii 



(24) 



The ratio of /hcii to /hi is related to r\ through r\ = 
(^Hc/^h)(/hoii//hi), where n Hc /n H « 1/13 (Y = 0.235) 
is the number ratio of helium to hydrogen atoms. Since Hen 
is a hydrogen-like species, the recombination coefficients 
scale similarly with temperature. Over the range 10,000 K 
to 20,000 K, QfHcin — 5.3qhii- Similarly, the photoionisa- 
tion rates also scale but in a more complicated manner de- 



pendent on the spectrum (Eq. [2}. The cross sections for 
Hi and Hen can be approximated by using Eq. [6] with J3 
and s the same for hydrogen-like species, aneii = cthi/4, 
and z/Hcn = 4i/hi (Table [2}. Taking the radiation field to 
have the form J„ oc v~ a , the integral in Eq. [2] gives for 
the ratio of photoionisation rates Thi/Thch = 2 2a+2 . Com- 
bining this result with Eq. [24] and aHein — 5.3«h 11 gives 
x] ~ 5.3/13 x 2 2a+2 . 

The values of r\ found in the simulations fluctuate about 
this estimate. For PL20, for example, the volume- averaged 
<V>voi = 11 at z = 3. However, the distribution is highly 
skewed; the mode and 68 percentile range is 3.31q i- The 
derived relation predicts 77 ~ 3.2 for a — 0.5, in good 
agreement with the mode. For MQ20, <rj> vo i — 14 and, 
like PL20, the distribution is highly skewed with the mode 
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Figure 14. The evolution of the distribution of r\ = »7.Hcii/ n Hi 
in the hybrid model HY20, for regions in which hydrogen was 
reionised by z = 6. Although the values of rj would be shifted by a 
change in the spectral shape of the ionising radiation background, 
the relative spread is invariant. The source spectrum makes a 
transition from a starburst spectrum to a power law between z = 
5 and z = 4. Because of the lack of He II- ionising photons, so that 
helium is not ionised to He III the distribution of r\ is much wider 
at z > 5 (not shown). 

3.8l u ' l 1 . The T] distribution is insensitive to redshift for the 
power-law and miniquasar models. In the case of the hybrid 
model, the evolution of the r/ distribution is complex after 
the power-law source turns on. For HY20, the evolution of 
the rj distribution is illustrated in Fig. [14] As the Hem I- 
front sweeps through the volume, r\ surges to > 1000. At 
the same time, the competing effect of the hardening spec- 
trum drives the most probable 77 values down. After the 
source spectrum has completed the transition to a power 
law, the most probably r\ value settles at the expected value 
of rj ~ 3.2 while the shadowed regions provide a high-77 tail. 

The fluctuations that produce the high-end tail of the 
77 distribution at z — 3 arise from the attenuation of the 
source radiation field at the Hen ionisation threshold. Re- 
call that we have corrected the data to remove /hi fluctu- 
ations, hence softer (larger values of a) local ionising fields 
are generated by attenuation. Variations in the local density 
have only a small role as they modify the local gas tem- 
perature which varies the approximation aHom — 5.3aHn- 
Figure [15] illustrates the interplay of effects. The bulk of the 
gas resides along the /Heii//m — 5.3 x 2 2a+2 locus in the 
/hoii — /hi plane. Attenuation leads to shadows with larger 
/hoii fractions, creating parallel loci of constant fm 11 //hi- 
The local overdensity sets the value of /hi, as we have set 
Thi to a constant value. Because of the correction to tihi, 
the effects of self-shielding are not seen in the 77 distribu- 
tion. Self-shielding would harden the spectrum, generating 
77 values below the bulk of the gas. This effect is seen in the 
uncorrected simulation data in the few regions in which it 
occurs. 

The magnitude of rj is sensitive to the input source spec- 
tra. For instance, choosing a — 2 instead of 0.5 for the 
power-law spectrum would increase 77 by a factor of 8 for 
the power-law model to values of 77 « 300. Boosting the ra- 
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Figure 15. Hell versus Hi fractions showing dependency on den- 
sity at z = 3 for model PL20, but with a uniform Hi photoionisa- 
tion rate (see text). The MQ20 and HY20 results are qualitatively 
similar while those for SB20 are uncorrelated owing to the paucity 
of helium-ionising photons. Only regions for which /hi < 0.1 by 
z = 6 are shown. 



tio of the contribution of the galaxy to the QSO spectrum in 
the hybrid model would achieve a similar effect. The relative 
spread in the fluctuations of 77, however, is independent of 
any overall shift in the amplitude or shape of the ionising 
background for regions in which both hydrogen and helium 
are ionised. The spread may be a useful diagnostic of the 
ionisation state of the IGM prior to complete Hen reioni- 
sation. Large fluctuations are found for 77 (or, equivalently, 
in /Heii//m). At z — 3, 77 is mostly constant, but there is 
small fraction with a large range (~ 2 dex). The extent of 
the range results both from the inhomogeneities in the ra- 
diation field and a wide spread in gas temperatures due to 
radiative transfer, particularly in the low density gas which 
gives rise to m ost of the Hen featu res. The spread is smaller 
than found bv lZheng et alj l|2004l ). who report a range of at 
least 2.5 dex in 77 at z < 3. This may indicate the presence 
of l ocal Hi- ionisi n g sour ces. 

IZheng et all l|2004) and iReimers etakl l|2005l ) report 
some regions with 77 > 1000, suggesting these may be re- 
gions for which Hen is still not ionised to Hem. Figure [TBI 
is a map of 77 for the hybrid model HY20 at 2 = 4 compared 
with a map of the /hoii fraction. Large 77 values are found to 
correspond to regions of high Hen fractions. The evolution 
of the 77 distribution for the hybrid model HY20 is illustrated 
in Fig. 1141 At z £ 4, when the power-law source is just tak- 
ing over from the starburst in magnitude, high excursions 
are found for 77, with values reaching up to 3000. These high 
77 regions are associated with incomplete Hen reionisation, 
where /hoii ;$ 1. Lowering the redshift at which the power- 
law source overtakes the starburst would lower the redshift 
at which thes e high excursions occ ur to z < 3. These results 
suggest that IReimers et alj ((2005) may have detected the 
Hen reionisation epoch. 

In their simulations, Ma selli fc Ferraral (|2005h find a 
higher value of 77 ~ 250, but a smaller spread of only a factor 
of two to three, not including the low values of 77 associated 
with self-shielding. As discussed earlier, our single I-front 
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Figure 16. Comparison of the Hell ionisation fraction (left panel) and the ratio of Hell to Hi number densities r) = »iHcii/ ra Hi (right 
panel), for the hybrid model HY20 at z = 4. (Regions for which hydrogen is not ionised by z = 6 are masked as blue.) A wide range in i) 
is found. The largest values (rj > 1000) correspond to regions where Hell is still incompletely reionised, with /h c ii ~ 1- The Hi fractions 
have been corrected to a uniform photoionisation rate (see text). 



simulations produce an over-abundance of shielded Hi gas 
at z < 6, which we have eliminated by correcting the hhi val- 
ues to a uniform Hi photoionisation rate. Their higher value 
of r? results partia l ly fro m the soft spectrum adopted from 
lHaardt fc Madaul (| 19961 ) for QSOs with a a = 1.5 source 
spectrum. 

We parametrize the /hi by the gas density, /hi = 
Q(p/p)' 3 , and similarly for /Heii- The best fit values are pro- 
vided in Tables [5] and [6] Cosmic variance leads to errors in a 
of about 10 per cent and /3 of about 2 per cent. Although the 
values of a may be readjusted by changing the magnitude 
of the ionisation radiation background, the values of j3 are 
invariant for a given source spectrum. The close agreement 
between the values of j3 for /hi and /hoii for the respective 
models reflects the near linear dependence between the Hi 
and Hen ionisation fractions, with a weak residual depen- 
dence on density (cf Fig. I15|l . 

The power-law dependence of ionisation fraction on 
density is expected if the shape of the local ionising spec- 
trum is constant and the density and temperature are re- 
lated through a power law. The equilibrium rate equation for 
hydrogen is nHiTHi = n e nHiiCtHii. Assuming almost com- 
plete ionisation, the ionisation fraction is, 



We recall that the photoionisation rate, Thi, is a function 
of the shape of the local ionising spectrum. Assuming it is 
constant and that ionisation is almost complete, the neutral 
fraction follows /hi oc pan u . Over the temperature range of 
interest, the recombination rate responds to the temperature 
as cchii oc T~ 69 . If the density and temperature are related 
polytropically, T oc p 7_1 , then /hi oc p 1-0 - 69 ^- 1 ) . From our 
power- law fits in Table [5] the inferred values for 7 at z — 3 
are 1.547± 0.001, 1.566 ±0.001, and 1.139±0.001 for PL20, 



MQ20, and HY20 respectively. For Hen the inferred values 
at z = 3 are 7 = 1.554 ± 0.001, 1.555 ± 0.001, and 1.222 ± 
0.001 PL20, MQ20, and HY20 respectively. The values of 
7 are nearly identical to those found by directly fitting the 
p-T relation in P31 



5 DISCUSSION AND CONCLUSIONS 

We have coupled a radiative transfer code based on a 
probabilistic photon transmission algorithm to a Particle- 
Mesh Af-body code in order to study the sensitivity of 
the post-ionisation temperature of the Intergalactic Medium 
on source spectrum. We performed multiple simulations 
with different spectra of ionising radiation: a power-law 
(oc i/ -0,5 ), miniquasar, starburst, and a time-varying hybrid 
spectrum that evolves from a starburst spectrum to a power 
law. 

The power-law and miniquasar spectra produce almost 
identical temperature distributions, owing to their similar 
shapes. A greater difference of temperatures is found be- 
tween the remaining models. The mass-weighted mean gas 
temperatures at z = 3 are 9000 K for the starburst source, 
15000 K for the power-law and miniquasar source, and 
17000 K for the hybrid models. A larger difference is found 
between the power-law/miniquasar and hybrid model tem- 
peratures from the polytropic fits to temperature and den- 
sity. A fit to the polytropic relation, T = To(p/p) 7_1 , gives 
T = 14600 ± 200, 14500 ± 160, 3515 ± 11, 18600 ± 100 and 
7 = 1.58 ± 0.01, 1.57 ± 0.01, 1.553 ± 0.002, 1.134 ± 0.005 
for the power-law, miniquasar, starburst and hybrid mod- 
els, respectively, at z = 3. The errors are formal fit errors. 
The highest temperatures are found in the hybrid model at 
the end of the transition from starburst to power-law spec- 
tra (2 ft* 4), at which time the temperatures span values up 



Reionisation scenarios and the temperature of the IGM 19 



Table 5. Fits to the relation /hi = a(p/p)@ for the gas for which /hi < 0.1 by z = 6. The errors are to the fits and not indicative of 
the distribution in the p — /h i plane. 

PL20 MQ20 HY20 



z 


a [xlO" 4 ] 


P 


a [xlO" 4 ] 


P 


a [xlO" 4 J 


P 


6 


5.093 ± 0.003 


0.847 ±0.001 


3.674 ±0.001 


0.8044 ± 0.006 


2.014 ±0.002 


0.803 ± 0.002 


5 


3.770 ± 0.002 


0.756 ± 0.001 


2.784 ± 0.001 


0.721 ± 0.001 


1.680 ± 0.001 


0.726 ±0.001 


4 


2.800 ± 0.002 


0.661 ±0.001 


2.083 ± 0.001 


0.643 ± 0.001 


1.8587 ± 0.001 


1.066 ±0.001 


3 


2.283 ±0.001 


0.6226 ± 0.0005 


1.7239 ± 0.0004 


0.6096 ± 0.0003 


2.022 ± 0.002 


0.904 ± 0.001 



Table 6. Fits to the relation /Hen = ot(p/p)^ for the gas for which /hi < 0.1 by z = 6. The errors are to the fits and not indicative of 
the distribution in the p — /hi plane. Note there is negligible Hen in the HY20 run prior to z = 5. 

PL20 MQ20 HY20 



2 


a [xl0 _/ ] 


P 


a [xlO -2 ] 


P 


a [xl0 _1! ] 


P 


6 


1.959 ± 0.002 


0.775 ± 0.001 


1.717 ± 0.002 


0.7508 ± 0.0009 






5 


1.531 ± 0.001 


0.695 ± 0.001 


1.3621 ± 0.0009 


0.6910 ± 0.0009 






4 


1.1989 ±0.0006 


0.6464 ± 0.0007 


1.0514 ±0.0004 


0.6380 ± 0.0005 


0.8657 ± 0.0008 


0.978 ± 0.001 


3 


1.0107 ± 0.0005 


0.6175 ± 0.0006 


0.8897 ± 0.0003 


0.6171 ± 0.0005 


0.9295 ± 0.0008 


0.847 ± 0.001 



to 40, 000 K , as required by measure ments of Civ and S iiv 
absorption systems in the Lya forest l|Rauch et al.| [l996). 

The He in I- front passing through Hn-dominated gas 
leads to a larger jump in temperature than the Hn I- front 
passing through Hem-dominated gas. Indeed, in the simula- 
tions with the power-law spectrum the temperature increase 
is only about 10 per cent when a trailing Hn I- front passes 
through gas in which helium is fully ionised. The difference 
is explained by the decrease in the mean molecular weight 
when Hi is ionised. Hence, a significant temperature change 
will not be produced around a hard source by the passage 
of a Hn I-front, in contrast to a soft source or a region in 
which hydrogen was fully ionised prior to helium. 

The post-ionisation temperature of the IGM may be 
used as a key observable for identifying the nature of the 
sources of reionisation. While moderate to high overdensity 
gas establishes an equilibrium temperature in which pho- 
toionisation heating balances atomic radiative cooling pro- 
cesses, the equilibrium time scale exceeds a Hubble time 
in lower density regions, those that give rise to optically 
thin Lya absorption systems. In contrast to optically thin 
reionisation models, we find a broad fanning out of the 
temperature-density relation for underdense regions, with 
temperatures exceeding 3 x 10 4 K for the hybrid model with 
late Hen reionisation. This may help substantially in recon- 
ciling the much larger Doppler parameters measured, which 
have a median value of about 30 km s _1 for the optically thin 
Lya systems at z fa 3, wit h those predicted b y simulations 
without radiative transfer (tMeiksin et al.lr200l| ). Even allow- 
ing for bulk motions to contribute as much as half to the line 
widths (in quadrature) , the large measured median Doppler 
parameters require a gas temperature of 2.8 x 10 4 K. The hy- 
brid model comes closest to meeting this requirement. Ulti- 
mately, incorporating hydrodynamics is necessary for defini- 
tive predictions of the detailed post-ionisation state of the 



IGM and a precise determination of the resulting statistical 
properties of the Lya forest. 

We do not address the effect of the epoch of reionisa- 
tion on the evolution of the gas temperature. Much earlier 
reionisation redshifts for the entire IGM will result in much 
cooler temperatures, primarily due to intense Compton cool- 
ing at high redshift. The effect of the reionisation epoch 
on the subsequent I GM t e mperature has been explo red by 
Ikeuchi fc Ostrikerl (Il986l ). iGiroux fc Shapiro! l| 19961 ). and 
Theuns et alj (|2002T ) without radiative transfer. A degener- 
acy exists between the epoch of reionisation and the nature 
of the sources on the temperature of the IGM. If the epoch 
of reionisation were determined through some other means, 
as for instance its detection in the radio by LOFAR or the 
MWA, the subsequent temperature of the IGM, particularly 
as revealed by the widths of optically thin Lya forest ab- 
sorbers, could then be used to determine the nature of the 
sources. 

Hardening of the spectrum due to passage through 
structures with high Hi column densities produces fluctua- 
tions in the /Heii//Hi ratio in the shadowed regions behind 
He in I- fronts prior to complete Hen rei onisation. A spread 
is indeed found in the data at z £j 3 (I Zheng et al 1 12004 
iReimers et afl l2005). The observed spread of about 2.5 dex, 
however, exceeds that found in our simulations at z — 3 by 
about a fact or of three. In parti cular, values of r\ > 1000 ar e 
reported by IZheng et alj i|2004h and IReimers et all l|2005l ). 
It may be that local sources are required to reproduce the 
wider range of observed values. Alternatively, at z ^ 4 val- 
ues up to r\ ~ 3000 are found for the hybrid model, when 
the power-law source begins to dominate over the starburst 
producing partially ionised Hen(3| Lowering the redshift at 



7 A similar effect is likely to occur for a power-law spectrum 
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which the power-law spectrum dominates over the starburst, 
thus delaying the Hen reionisation epoch, would lower the 
redshift at which large fluctuatio ns in rj are produced. It 
may be that lReimers et al.l (|2005l ) have detected the epoch 
of Hen reionisation. 

Not explored in detail by this paper is the pre-ionisation 
state of the gas which is relevant to future 21 cm detections. 
To correctly model the gas will require the incorporation of 
the production of secondary electrons and the diffuse radia- 
tion field. The production of the secondary electrons further 
cools the gas due both to overcoming the ionisation potential 
and to radiative losses fro m the conse quent enhancement of 
Lya collisional excitation (|Shul]Hl979T l. A copious number of 
Lya photons could be produced in any region where helium 
is ionised prior to hydrogen. The photons may have impor- 
tant implications for the detection of the 21cm signature of 
the IGM, as they can decoup le the spin temperature from 
the CMB (jMadau et aflTjjj ). The Lya photons will also be 
an additional source of pre-heating through recoils, although 
the amount i s likely to be s mall before the scattering reaches 
equilibrium (|Me iksin 200^). We are currently including both 
hydrodynamics and the extra radiative physics in more re- 
fined models. The addition of hydrodynamics will also elim- 
inate any limitations owing to the absence of advection. 
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